function grad = ComputeGradient(Function, rs)    

%% Check 
    eps = 1.0;      
    %% first we find the smallest eps we can handle
    one = 1.0;
    two = 2.0;
    test = 0.0;
    
    nDim = size(rs);
    nDim = nDim(1);
    
    while (1)
        eps = eps/two;
        test = eps * one + one;
        if (test > 1.0)
            continue;
        else 
            break;
        end
    end    
    eps = eps*two*two;
    
    %% now let's go for the finite different        
    grad = zeros(size(rs));
    fx = feval( Function, rs);       
    for i=1 : nDim,        
      grad(i) = GradComponent(Function, rs, i, eps, fx);
    end   
end 